home *** CD-ROM | disk | FTP | other *** search
/ GFX Sensations 1 / Graphic Sensations - Volume 1.iso / tools / amiga / 3d_tools / irit40s.lha / Irit / cagd_lib / cagdoslo.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-30  |  7.7 KB  |  208 lines
  1. /******************************************************************************
  2. * CagdOslo.c - an implementation of the oslo algorithm (1).              *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Aug. 92.                          *
  5. ******************************************************************************/
  6.  
  7. #include <cagd_loc.h>
  8.  
  9. #define KNOT_INFINITY 1e20
  10.  
  11. #ifdef DEBUG
  12. static void PrintAlphaMat(BspKnotAlphaCoeffType *A);
  13. #endif /* DEBUG */
  14.  
  15. /******************************************************************************
  16. * Computes the values of the alpha coefficients, Ai,k(j) of order k:          *
  17. *                                                                             *
  18. *         n              n    m                   m    n              *
  19. *         _              _    _                   _    _              *
  20. *  C(t) = \ Pi Bi,k(t) = \ Pi \ Ai,k(j) Nj,k(t) = \  ( \ Pi Ai,k(j) ) Nj,k(t) *
  21. *         /              /    /                   /    /              *
  22. *         -              -    -                   -    -              *
  23. *        i=0            i=0  j=0                 j=0  i=0              *
  24. *                                                                             *
  25. * Let T be the original knot vector and let t be the refined one, i.e. T is   *
  26. * a subset of t.                                   *
  27. * The Ai,k(j) are computed from the following recursive    definition:           *
  28. *                                                                             *
  29. *             1, T(i) <= t(i) < T(i+1)                                        *
  30. *           /                                                                 *
  31. * Ai,1(j) =                                                                   *
  32. *           \                                                                 *
  33. *             0, otherwise.                                                   *
  34. *                                                                             *
  35. *                                                                             *
  36. *           T(j+k-1) - T(i)             T(i+k) - T(j+k-1)                     *
  37. * Ai,k(j) = --------------- Ai,k-1(j) + --------------- Ai+1,k-1(j)           *
  38. *           T(i+k-1) - T(i)             T(i+k) - T(i+1)                       *
  39. *                                                                             *
  40. * LengthKVT + k is the length of KVT and similarly LengthKVt + k is the       *
  41. * length of KVt. In other words, LengthKVT and LengthKVt are the control      *
  42. * points len...                                      *
  43. *                                                                             *
  44. * The output matrix has LengthKVT rows and LengthKVt columns (#cols > #rows)  *
  45. * ColIndex/Length hold LengthKVt pairs of first non zero scalar and length of *
  46. * non zero values in that column, so not all LengthKVT scalars are blended.   *
  47. ******************************************************************************/
  48. BspKnotAlphaCoeffType *BspKnotEvalAlphaCoef(int k, CagdRType *KVT,
  49.                 int LengthKVT, CagdRType *KVt, int LengthKVt)
  50. {
  51.     int i, j, o;
  52.     CagdRType *m, **Rows;
  53.     BspKnotAlphaCoeffType
  54.     *A = (BspKnotAlphaCoeffType *)
  55.                      IritMalloc(sizeof(BspKnotAlphaCoeffType));
  56.  
  57.     A -> Order = k;
  58.     A -> Length = LengthKVT;
  59.     A -> RefLength = LengthKVt;
  60.     A -> Matrix = (CagdRType *) IritMalloc(sizeof(CagdRType) * (LengthKVT + 1)
  61.                                  * LengthKVt);
  62.     Rows = A -> Rows = (CagdRType **) IritMalloc(sizeof(CagdRType *) *
  63.                                   (LengthKVT + 1));
  64.     A -> ColIndex = (int *) IritMalloc(sizeof(int) * LengthKVt);
  65.     A -> ColLength = (int *) IritMalloc(sizeof(int) * LengthKVt);
  66.  
  67.     /* Update the row pointers to point onto the matrix rows. */
  68.     for (i = 0, j = 0; i <= LengthKVT; i++, j += LengthKVt)
  69.     Rows[i] = &A -> Matrix[j];
  70.  
  71.     /* Initialize the matrix with according to order 1: */
  72.     for (m = A -> Matrix, i = (LengthKVT + 1) * LengthKVt; i > 0; i--)
  73.     *m++ = 0.0;
  74.     for (i = 0; i < LengthKVT; i++) {
  75.     CagdRType
  76.         *RowI = Rows[i],
  77.         KVTI = KVT[i],
  78.         KVTI1 = KVT[i + 1];
  79.  
  80.     for (j = 0; j < LengthKVt; j++, RowI++) {
  81.         if (KVTI <= KVt[j] && KVt[j] < KVTI1)
  82.         *RowI = 1.0;
  83.     }
  84.     }
  85.  
  86.     /* Iterate upto order k: */
  87.     for (o = 2; o <= k; o++) {
  88.     for (i = 0; i < LengthKVT; i++) {
  89.         CagdRType
  90.         *RowI = Rows[i],
  91.         *RowI1 = Rows[i + 1],
  92.         KVTI = KVT[i],
  93.         KVTIO = KVT[i + o],
  94.         t1 = KVT[i + o - 1] - KVTI,
  95.         t2 = KVTIO - KVT[i + 1];
  96.  
  97.         /* If ti == 0, the whole term should be ignored. */
  98.         if (t1 < EPSILON)
  99.         t1 = KNOT_INFINITY;
  100.         if (t2 < EPSILON)
  101.         t2 = KNOT_INFINITY;
  102.  
  103.         for (j = 0; j < LengthKVt; j++, RowI++, RowI1++) {
  104.         int jo = j + o;
  105.  
  106.         *RowI = *RowI * (KVt[jo - 1] - KVTI) / t1 +
  107.             *RowI1 * (KVTIO - KVt[jo - 1]) / t2;
  108.         }
  109.     }
  110.     }
  111.  
  112.     /* Update the row non zero indices. */
  113.     for (i = LengthKVt - 1; i >= 0; i--) {
  114.     for (j = 0; j < LengthKVT && APX_EQ(Rows[j][i], 0.0); j++);
  115.     A -> ColIndex[i] = j;
  116.     for (j = LengthKVT - 1; j >= 0 && APX_EQ(Rows[j][i], 0.0); j--);
  117.     if ((A -> ColLength[i] = j + 1 - A -> ColIndex[i]) <= 0)
  118.         FATAL_ERROR(CAGD_ERR_DEGEN_ALPHA);
  119.     }
  120.  
  121.     return A;
  122. }
  123.  
  124. /******************************************************************************
  125. * Same as EvalAlphaCoef but the new knot set NewKV is merged with KVT to      *
  126. * form the new knot vector KVt.                              *
  127. ******************************************************************************/
  128. void BspKnotFreeAlphaCoef(BspKnotAlphaCoeffType *A)
  129. {
  130.     IritFree((VoidPtr) A -> Matrix);
  131.     IritFree((VoidPtr) A -> Rows);
  132.     IritFree((VoidPtr) A -> ColIndex);
  133.     IritFree((VoidPtr) A -> ColLength);
  134.     IritFree((VoidPtr) A);
  135. }
  136.  
  137. /******************************************************************************
  138. * Same as EvalAlphaCoef but the new knot set NewKV is merged with KVT to      *
  139. * form the new knot vector KVt.                              *
  140. ******************************************************************************/
  141. BspKnotAlphaCoeffType *BspKnotEvalAlphaCoefMerge(int k, CagdRType *KVT,
  142.                    int LengthKVT, CagdRType *NewKV, int LengthNewKV)
  143. {
  144.     int LengthKVt;
  145.     BspKnotAlphaCoeffType *A;
  146.  
  147.     if (NewKV == NULL || LengthNewKV == 0) {
  148.     A = BspKnotEvalAlphaCoef(k, KVT, LengthKVT, KVT, LengthKVT);
  149.     }
  150.     else {
  151.     CagdRType
  152.         *KVt = BspKnotMergeTwo(KVT, LengthKVT + k,
  153.                    NewKV, LengthNewKV, 0, &LengthKVt);
  154.  
  155.     A = BspKnotEvalAlphaCoef(k, KVT, LengthKVT, KVt, LengthKVt - k);
  156.  
  157.     IritFree(KVt);
  158.     }
  159.  
  160.     return A;
  161. }
  162.  
  163. /******************************************************************************
  164. * Prepare a refinement vector for the given knot vector domain with n          *
  165. * inserted knots equally spaced.                          *
  166. ******************************************************************************/
  167. CagdRType *BspKnotPrepEquallySpaced(int n, CagdRType Tmin, CagdRType Tmax)
  168. {
  169.     int i;
  170.     CagdRType dt, t, *RefKV;
  171.  
  172.     if (n <= 0) {
  173.     FATAL_ERROR(CAGD_ERR_WRONG_INDEX);
  174.     return NULL;
  175.     }
  176.  
  177.     dt = (Tmax - Tmin) / (n + 1),
  178.     t = Tmin + dt,
  179.     RefKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * n);
  180.  
  181.     for (i = 0; i < n; i++, t += dt) {
  182.     RefKV[i] = t;
  183.     }
  184.     return RefKV;
  185. }
  186.  
  187. #ifdef DEBUG
  188. /******************************************************************************
  189. * Prints the content of the alpha matrix.                      *
  190. ******************************************************************************/
  191. static void PrintAlphaMat(BspKnotAlphaCoeffType *A)
  192. {
  193.     int i, j;
  194.  
  195.     fprintf(stderr, "Order = %d, Length = %d\n", A -> Order, A -> Length);
  196.     for (i = 0; i < A -> Length; i++) {
  197.     for (j = 0; j < A -> RefLength; j++)
  198.         fprintf(stderr, " %9.5lg", A -> Rows[i][j]);
  199.     fprintf(stderr, "\n");
  200.     }
  201.  
  202.     fprintf(stderr, "    ");
  203.     for (j = 0; j < A -> RefLength; j++)
  204.         fprintf(stderr, " %3d %3d |", A -> ColIndex[j], A -> ColLength[j]);
  205.     fprintf(stderr, "\n");
  206. }
  207. #endif /* DEBUG */
  208.